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The process of equilibration in theory is investigated for a homogeneous system in 3+1 
dimensions and a variety of out-of-equilibrium initial conditions, both in the symmetric and 
broken phase, by means of the 2PI effective action. Two d>-derivable approximations includ¬ 
ing scattering effects are used: the two-loop and the “basketball”, the latter corresponding 
to the truncation of the 2PI effective action at 0{X^). The approach to equilibrium, as well 
as the kinetic and chemical equilibration is investigated. 


I. INTRODUCTION 

The approach to equilibrium is an important aspect of non-equilibrium dynamics. In the context 
of particle physics, a large part of the interest derives from results of heavy-ion collision experiments 
with the RHIC at Brookhaven. The hydrodynamic description of the experiments suggests that 
there is early thermalization [1], but a short thermalization time seems to contradict traditional 
perturbative estimates [2, 3]. This puzzle has been analyzed in terms of prethermalization [4], and 
led to further study of the microscopic dynamical processes responsible for the equilibration of the 
quark-gluon plasma [5? ? , 6]. Understanding the dynamical processes leading to equilibration in 
theories with simpler interactions may also shed some light on this issue. We focus in this paper 
on the case of scalar theory. 

An adequate method to study out-of-equilibrium dynamics from first principles is the closed- 
time-path formalism [7-10]. This scheme leads to causal equations of motion for the various 
correlation functions that describe their time evolution. The initial conditions are specified by 
a density matrix, which can be far from equilibrium. Most applications of this method have 
focused on the study of the equations of motion for the 1- and 2-point functions, known as the 
Kadanoff-Baym equations [11]. Close enough to equilibrium, where kinetic theory is applicable, 
the Kadanoff-Baym equations have been used extensively, mostly in connection with the study of 
transport phenomena and the derivation of effective Boltzmann equations (see, for instance [12-17] 
and references therein). 

Far from equilibrium, kinetic theory is no longer valid, and simple perturbation theory ap¬ 
proaches fail to work due to the appearance of secular terms (see for instance [18]) and/or pinch 
singularities [19]. These problems are usually absent if one makes use of a self-consistent method, 
such as the Hartree approximation. Unfortunately, real-time Hartree descriptions do not include 
sufficient scattering between the field modes, and thus fail to describe the approach to equilibrium. 
They are also not “universal”, in the sense that the memory from the initial configuration is not 
completely lost [20]. An infinite number of conserved charges appear that prevent the system from 
reaching a universal equilibrium state, independently of the initial conditions. However, a Hartree 
ensemble approximation has been formulated to give an improved description of the early approach 
to equilibrium [21, 22]. 

When the particle occupation numbers are large, another useful method far from equilibrium 
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is the classical approximation. Interesting situations where this occurs include preheating after 
cosmological inflation due to parametric resonance [23-26] or spinodal decomposition, [23, 27-31] 
as well as the early stages of a heavy-ion collision [32-35], where the gluon occupation numbers are 
as large as ~ l/a^, up to a saturation scale [2, 36]. The classical approximation is not good for 
describing quantum equilibration, since the system does not move towards the quantum, but to 
the classical equilibrium state. Nevertheless, the classical theory has been used to shed some light 
on the dynamics of equilibration and relaxation [30, 37-40], as well as a testground for comparison 
with various other approximation schemes [20, 31, 41, 42]. 

A powerful scheme that takes into account both scattering and quantum effects is the two- 
particle irreducible (2PI) effective action [14, 43, 44]. The 2PI effective action furnishes a complete 
representation of the theory in terms of the dressed 1- and 2- point functions. The exact equations 
of motion describing the time evolution of these correlation functions are obtained by a variational 
principle on the 2PI effective action functional. Various approximations to the equations of motion 
can be obtained if one applies the variational method to a truncated version of the action. By 
construction, these are self-consistent and thus free of secular problems. The approximation can 
be improved, in principle, by truncating the 2PI effective action at higher order in some expansion 
parameter. 

The main advantages of the 2PI effective action approach stem from the fact that the approx¬ 
imations are performed on the level of a functional. For that reason the approximations have 
also been called Functional-derivable, or <h-derivable. The 2PI effective action functional (and any 
truncation thereof) is, by construction, invariant under global transformations of the 1- and 2-point 
functions. The variational procedure on any truncation guarantees that the global symmetries are 
still preserved by the equations of motion. Their associated Noether currents are thus conserved. 
In particular, this implies that the derived equations of motion conserve energy, as well as global 
charges [45-47]. This is a very important feature when studying out-of-equilibrium processes, 
where most other quantities evolve in complicated ways. The <h-derivable approximations to the 
2PI effective action constitute thus a very convenient method for studying equilibration. 

In recent years, approximations based on the 2PI effective action have been applied succesfully 
to the study of non-equilibrium real-time dynamics. In the context of scalar theories, studies of 
equilibration have been carried out in the 3-loop 4>-derivable approximation for theory, in the 
symmetric phase, both in 1-|-1 [44, 48] and 2-|-l dimensions [49]. In the broken phase, the case of 
1-|-1 dimensions has also been discussed in [50, 51]. Similar studies of thermalization have been 
performed in the 0{N) model in 1-|-1 dimensions, both at next-to-leading (NLO) order in a \/N 
expansion [52, 53], and in the bare vertex approximation [50, 51, 54]. All these analyses, which 
include scattering, show that the system indeed equilibrates, with the equilibrium state independent 
of the initial conditions. Comparing with the loop expansion, the 1/N expansion has the advantage 
that it is applicable in situations where large particle numbers are generated. This has allowed 
the study of interesting phenomena, such as parametric resonance [55] or spinodal decomposition 
during a phase transition [31]. The studies in [55] and [31] were done for the 0{N) model at NLO, 
in 3-|-l dimensions. Methods based on the 2PI effective action have also been applied to theories 
with fermions, in 3-|-l dimensions [4, 56]. The extension of the 2PI effective action methods for 
gauge theories, however, is not straightforward due to a residual dependence on the choice of gauge 
condition [57? -59]. 

In this paper, we use the loop-expansion of the 2PI effective action to study the approach to 
equilibrium for the case of a real scalar theory in 3-1-1 dimensions, both in the symmetric and 
broken phase, complementing in this manner the investigations in [44, 49, 51]. 
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II. 2PI LOOP EXPANSION OF THEORY 


For ip‘^ theory we write the action as 



]^d^ip{x)d^^ip{x) - ]Wip{xf 



( 1 ) 


The subscript C indicates that the integrations are performed along the real-time Schwinger- 
Keldysh contour, running from an initial time to to time t along C+ and going back to to along 
C- (see Fig. 1). The formulation of the theory along the real-time contour C is appropiate for 
studying non-equilibrium problems [10, 13, 60]. 

The system can be in two distinct phases: the symmetric phase (the vacuum field expectation 
value u is u = 0) which occurs for > 0, and the broken phase {v ^ 0) for rr? < 0. At tree level, 
the vacuum expectation value in the broken phase is given by utree = Y^6|m^|/A. 



j ^ Schwinger-Keldysh contour 

c. C 



^ J t 

0 Q t 


Figure 1: Schwinger-Keldysh contour. 


The complete information about the theory can be written in terms of the 2PI effective action, 
which depends explicitly on the full connected 1- and 2-point functions (/>(x) = {^{x)) and G{x, y) = 
{Tcip{x) (p{y)) — (j){x)(j){y). For scalar theory, the 2PI effective action functional can be written 
as [43] 


with 


mG] = S[^]-'-TrlnG+^-Tr 


(Go'-G-^)-G +mG], 


*Go^(x,y) 


64>{x)5(p{y) 




^c{x,y). 


The contour delta function 6c{x,y) is given by 


( 2 ) 

( 3 ) 


Sc{x,y) 


1 if X = y and x, y G C+, 

< —1 if X = y and X, y G C_, 
0 otherwise. 

\ 


( 4 ) 


The functional comprises the sum of the closed two-particle-irreducible (2PI) skeleton diagrams. 
Up to three loops it is given by 
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Truncation 

Order 

G] 

Hartree approximation 

0{\) 


Two-loop approximation 

2 loops 

1 0 -L 1 ,, ,, 

8 Q + 12 " “ 

“Basketball” approximation 

0[X^) 

8 0 + 12 ’‘O’* + 48 


Table I: Truncations of the 2PI effective action. 


In this manner, the functional ^[(j),G] is 

^[0,^] = -^ J G{x,xf+i^ j (fx J (l){x)G{x,yf(l){y)+i^ j d'^x J d'^y G{x,yf+... 

( 

The 2PI effective action r[(/i, G] provides an exact representation of the full theory. Considering 
only a finite number of terms in the series of diagrams in leads to a truncated action, from which 
approximate “physical” 1- and 2-point functions can be obtained by a variational procedure. As a 
result of this, a resummation of effects from higher orders in perturbation theory is performed. 

In this paper we shall investigate the truncations of the 2PI effective action up to three loops, 
in particular up to O(A^). The various truncations considered and their corresponding truncated 
functionals 4>tr are displayed in table I. The organization of the truncations discussed is based on 
the superficial counting of loops and/or vertices in the diagrams of <1>, i.e. no assumption is taken 
on the coupling constant dependence of cj) or G. 

In our analysis using the three-loop “basketball approximation” we have neglected the other 
three-loop diagrams 

These diagrams are respectively of superficial order 0{}?(lP‘) and 0(A^i/)^). In the symmetric phase, 
where (/> ~ 0, they can be safely neglected. In the broken phase, however, (/> ~ Utree ~ 0(A and 
thus both diagrams become O(A^). In this situation it is not clear whether these contributions can 
be ignored. Due to the difficulty in treating the above diagrams numerically, we decided to neglect 
them in our analysis. Part of the first diagram in (8) can be recovered at NLO in a 1/A^-expansion 
[53]. 


III. EQUATIONS OF MOTION 


In the formulation on the real-time contour C, a <h-derivable approximation to the 2PI effec¬ 
tive action T leads to equations of motion for the 1- and 2-point functions. Indeed, solving the 
stationarity conditions 


smG\ 

5(j) 


= 0 , 


leads to the equation for the mean field 


5G 


= 0 , 


SS[^] 

6(l){x) 


+ ^>^G{x,x)(j){x) = 


smG] 

6(j){x) 


(9) 


( 10 ) 
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and for the 2-point function 


^c{x,y)= / (fzGQ\x,z)G{z,y)+i / <fzT.{x,z)G{z,y). 

Jc Jc 

The self-energy S(x,y) is given, in terms of the functional by^ 

mcp,G] 


^{x,y) = -2 


6G{y,x)' 


( 11 ) 


( 12 ) 


To the order considered here, the self-energy S is determined from the truncated functional 
To 0(A2) one finds 


mG]=i 




(13) 


For the case of the Hartree approximation (see Table I), only the first diagram in (13) (the “leaf” 
diagram) enters in S. For the case of the two-loop and “basketball” approximations respectively, the 
second and third diagrams in (13) (the “eye” and the “sunset”) have to be taken into account. For 
the study of nonequilibrium dynamics these are important diagrams as they account for scattering 
and hence can lead to equilibration. 

The self-energy can be split up into a local and a nonlocal part, 

S(x, y) = S'(x)(ic(a:, y) + S°^(x, y), (14) 


with 

'^^{x) = ^G{x,x), (15) 

y) = -i^(p{x)G{x, yf(t){y) - i^G{x, yf. (16) 

2 D 

The quantities entering in the equations of motion (10) and (11) are defined in the real-time con¬ 
tour C. For the non-local quantities, such as G{x,y) and T,'^^{x,y), this implies the appearance 
of several components, corresponding to the various positions of the time indices along the con¬ 
tour. For G{x,y), the various contour components are written in a compact manner by using 
the decomposition in terms of the correlators G^{x,y) = {(p{x)(p{y)) and G^{x,y) = {ip{y)ip{x)), 
namely 


G{x,y) = 0c(a:o - yQ)G'^{x,y) + 0c(yo - a^o)G<(x,y). (17) 

The 0-functions used here are defined along the contour C. A similar decomposition can be written 
for the self-energy S“*(x,y), i.e. 

S“'(x, y) = 0c(xo - yQ)T?'{x, y) + Qc{yo - xo)S<(x, y). (18) 

From (17) we see that the dynamics of the propagator is entirely described by the two complex 
functions G^ and G^. For the real scalar theory under consideration, these functions satisfy the 
property [G^{x,y)]* = G^{x,y), which leaves only one independent complex function describing 


^ Our convention for the self-energy E is that it appears, formally, as a positive contribution to the mass. In 
particnlar, it is given in terms of the self-energy Es used in [52] by E = iEs. 








the propagator dynamics. This can be parametrized in terms of two real functions F and p 
according to 


G>{x,y) = F{x,y)-^p{x,y), (19) 

G^{x,y) = F{x,y)+^-p{x,y). (20) 

The functions F and p correspond to the correlators 

F{x:y) = ^ [G>{x,y) + G<{x,y)] = ]^{{‘p{x),ip{y)]), (21) 

p{x,y) = iG>{x,y) -iG^{x,y) = i {[ip{x),ip{y)]). (22) 


The correlators F{x, y) and p{x^ y) contain, respectively, statistical and spectral information about 
the system. They satisfy the symmetry properties 


F{x,y) = F{y,x), 
p{x,y) = -p{y,x), 

which make them very useful for numerical implementation [48]. 

For the self-energy we introduce, in a similar fashion, the quantities 


(23) 

(24) 


S^(x,y) = ^ [S>(x,y) + S<(x,y)] , 


Y.P{x,y) = Y.<{x,y) - S>(x,y), 


(25) 

(26) 


which satisfy similar properties as their propagator counterparts. For the case of the non-local 
part of the self-energy given by (16), these become 


S^(x,y) = ^<j){x)<j){y) 


F‘^{x,y) - 


p^{x,y) 


A2 


+ ^F(x,y) 1^2 

6 


F‘^{x,y) - 


3p‘^{x,y) 


T,P{x,y) = \^(t){x)4){y)[F{x,y)p{x,y)] + —p{x,y) 


3F^{x,y) - 


p^ix,y) 


(27) 

(28) 


In the study presented here, we shall focus on the dynamics of the statistical and spectral cor¬ 
relators F and p. Their equations of motion are determined from (11) by using the decompositions 
(17) and (18), as well as the definitions (21) and (22). For the case xq > yo, one finds 

f-xo f ryo r 

[dl +M^{x)] F{x,y) = J dzo J d^zT.P{x,z)F{z,y) - J dzo J d^z {x, z)p{y, z), (29) 


rxo 


[dl +M^{x)] p{x,y) = / dzo / d^z E^^ix, z)p{z,y), 


with 


'yo 


M^(x) = 171^ + ^(/>(x)^ -|- S^(x) = -|- ^(pix)"^ + ^F{x, x). 


(30) 


(31) 


With the same considerations as for the 2-point functions, the equation of motion of the mean field 
4>{x) is found from (10) to be 
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where Yl^(x,z) is the p-component of the “sunset” self-energy diagram, given by 

=-yp(x, 2 :) 3F{x, zf - . (33) 

This contribntion derives from inclnding the second of the 2PI diagrams in (see Eq. 5). Therefore 
it is present in both the two-loop and “basketball” approximations. The tilde in S is written to 
avoid any confusion with the self-energy S entering in the equations of motion for the propagator. In 
that case, the “sunset” diagram enters in the self-energy S only in the “basketball” approximation. 

With the self-energies and given respectively by (27),(28) and (33), equations (29-32) 

constitnte a set of closed coupled evolution eqnations for the correlators F and p and the mean 
field (f). These equations are explicitly cansal, i.e. the evolution of F, p and 4* is determined by 
the values of the correlators and mean fields at previons times. The driving terms in the RHS 
of those equations consist of nonlocal “memory” integrals that contain the information about the 
earlier stages of the evolution. By specifying a complete set of initial conditions for F, p and (/>, the 
equations of motion (29-32) constitute an initial value problem. We perform a numerical analysis 
of the equations of motion in the next section. 

We finish this section with the calculation of the energy density corresponding to the truncations 
of the 2PI effective action. The energy density is determined from the energy-momentum tensor 
component It takes the form (see appendix A) 

r°°(x, t) = ]- dtdt' + 5x9x' + (F(x, t;x', t') + (/>(x, t)(/>(x, t')) 

2 L J x=x' 

t=t' 

1 A \ O 

+ + ^A(/>(x, t)2F(x, t; X, t) - (34) 

Here ^( 2 ;) is an auxiliary scale factor introduced in the conpling constant as A ^ C(3j)A. For a 
given truncation, the energy density is obtained from (34) by the snbstitution ^ <I>tr. In the 
“basketball” approximation, for instance, the energy density becomes 

r°°(x,t) = ^ [at(/>(x,t)]^ ^ [5x<()(x,t) • 9x<(>(x,t)] -h + ^5x • 5yF(x,t;x,t)|^^y 

+ ]^rn^ t) + F(x, t; x, t)] -h ^A(()(x, + |AF(x, t; x, t)(()(x, tf + ^AF(x, t; x, tf 

+ y/ dzo J d^z - 3p{x,t;z, zo)F{x,t;z, zof 4>{z,zo) 

y / J - F{:>c,t;z,zof F{x,t-,z, zo)p{x,t;z, zo). (35) 

It follows from translational invariance [45, 46] (see also appendix A), that the energy density (35) 
is exactly conserved in the evolntion. 

IV. NUMERICAL ANALYSIS AND RENORMALIZATION 

We study the non-equilibrium evolution of the correlators F and p and the mean field (j) by 
solving numerically the equations of motion (29-32), both in the symmetric and the broken phase. 



A. Numerical implementation 


We shall consider the system to be discretized on a space-time lattice with a finite spatial volume 
and spatially periodic boundary conditions. The action of theory on a space-time lattice is 


‘S'latl'/?] = a^at'^ 

x,i 




(36) 


The lattice spacings a and at correspond to the spatial and time directions, respectively. The 
derivatives stand for forward finite differences, e.g. t) = (l/at) [</?(x, t -|- at) — The 

spatial lattice volume is given in terms of the number of lattice sites N as V = = {Na)^. In 

the following we shall use lattice units (a = 1) and write dt = at/a for the dimensionless time-step. 
The mass mo and coupling Aq are hare parameters to be determined below. The lattice version of 
the squared spatial momentum is given by 

3 27r7?/' ./V TV 

^ (2 - 2 cos kj), with kj = ni = + 1,...,— (for iV even). (37) 

i=l 


Plotting data as a function of corrects for a large part of the lattice artifacts. 

The lattice provides a cutoff and regularizes the ultraviolet divergent terms in the continuum 
limit, which are to be dealt with by renormalization. The continuum renormalization of <I>-derivable 
approximations has been studied in detail in [61-65]. For our purpose it is enough to use an 
approximate renormalization that ensures that the relevant length scales in our simulations are 
larger than the lattice spacing a. This is achieved by simply choosing the bare paremeters mo and 
Aq according to the one-loop formulas that relate them to the renormalized parameters. The bare 
mass mo is given in terms of the renormalized mass m by (see also [63, 66, 67]) 

mp = m^ — 5m?, (38) 


with the mass counterterm 


5rr? 


2a2 


/i(am) 


A2u2 

2a? 


hiam). 


(39) 


The Ii and I 2 in (39) are dimensionless integrals coming respectively from the one-loop “leaf” and 
“eye” diagrams in the self-energy S at zero temperature. On the lattice, for continuous time and 
in the infinite volume limit (N 00 ), they are given by 


/i(am) = i 




dkn 


d^k 


-n/a (2vr)3 J 2Tt kl-a -rr? -ie J (27r)3 


rf,- r" d^k f dko ( 


1 


d^k 


lat 

1 


'-^/a(2vr)3 7 2tt \kl-a-^kl^-m?-iej J_^ {a^rr? + kl,)' 

The bare coupling Ao can be determined from the one-loop expression (see [63, 66, 67]) 


(40) 

:■ (41) 


1 

Ao 


- - hiam). 


(42) 


The renormalization conditions that define the renormalized mass and coupling via Eqns. (38) 
and (42) are such that they correspond to the values of the two- and four-point vertex functions 
at vanishing external momenta. For the values of the couplings (A = 1,6) and lattice spacing 
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(0.5 < am < 1) that we use in our simulations, the difference between Aq and A is less than 10%. 
In practice, we simply choose Aq as if it were the renormalized coupling. 

For the renormalization of the mass we use (38), which for the case of a spatial lattice is 
then given by 


2 I 2 

ruQ = ±m 


^0 _1_ i_ _1_ 

4a2iV3 ^ VaW + Zci;, ^ V(aW + fey' 


(43) 


In practice, we conveniently choose a value for the renormalized mass (such that am < 1), which 
determines via Eq. (43) the bare mass that enters in the equations of motion for the mean field 
and propagator. In our simulations we used am = 0.7 and N = 16. Given the input parameters 
mo and Aq, the output physics is of course not known precisely, it is determined by the <I>-derived 
equations of motion. 


B. Initial conditions 

We specialize to a spatially homogeneous situation. In this case, F{x,y) = F{t,t',x — y), 
p{x,y) = p(t,t'— y) so we can perform a Fourier transformation and study the propagator 
modes Pk{t,t') and Fk(t,t'). In addition, the mean field depends only on time. To specify the 
time evolution, the equations of motion (29-32) must be supplemented with initial conditions at 
t = t' = 0. These are given by the values and derivatives of p, F and (p at initial time. The initial 
conditions for p follow from it being the expectation value of the commutator of two fields, which 
implies 


Pk{t,t) = 0 , dtpk{t,t')\^^^, = I, 


Imposing the condition (44) at t = t' = 0, it is preserved by the equations of motion. 
For the statistical correlator F, we choose initial conditions of the form 


({7rk(t),7?-k(t')})lt=t'=o = 9iTk(t,tO|i=t/=o 
({7rk(t)7r_k(t0})lt=t'=o = 9t9t/Fk(t,tO|i=t/=o 


1 

tUk 




= 0 , 


= Wk 


nk + i 


(44) 


(45) 

(46) 

(47) 


where 7rk(t) = dt(p\i{t) are the conjugate field momenta, Uk is some distribution function and tUk = 
%mfjj -|- , with min to be specified shortly. An initial condition of this form can be represented 

by a gaussian density matrix. We will use the following cases for the distribution function nk: 


a) Thermal: The distribution function nk corresponds to a Bose-Einstein, at some initial tem¬ 
perature Tin, 

^ e%k/Tn) _ 1 ■ 

This also includes the “vacuum” initial condition of Tjn = 0 (which is of course only an 
approximation to the vacuum state in the interacting theory). The input mass for all Tin is 
the renormalized min = m in the symmetric phase, and mjn = \/2m in the broken phase. 
We expect these to be close to the zero-temperature particle masses, respectively in these 
phases. 
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b) Top-hat: In this case, only modes with momenta within a range < k ^^^ are 

occupied. The distribution function can be parametrized as 

nk = r? 0(kLx - k2)0(k2 - (49) 


where rj represents the occupancy of the excited modes. The input mass is again given by 
ruin = m (symmetric) and m\n = (broken). 

The mean held is initialized at (/> = 0 (“symmetric phase”) and (p = utree; (the zero temperature 
“broken phase”). Below, we will also allow the mean held to be slightly displaced from these two, 
in order to study relaxation in a thermal background. 


C. Observables 


As the system evolves in time, we expect the scattering processes to lead to equilibration. The 
occupation numbers of the momentum modes are expected to gradually approach a Bose-Einstein 
distribution, provided the coupling is not too strong. The statistical information about the evolving 
system can be extracted from the equal-time correlation function F\^{t, t). We can use Tk to dehne 
a quasiparticle distribution function and frequencies as [21, 22, 30, 38, 48] 


nUi) + 2 ^ F^{t,t), 

^p(t) — 


ldtdt'Fk{t,T] 


t=t' 


Fk(t,t) 


(50) 

(51) 


The correction Ck diminishes errors associated with the time discretization on the lattice. It is 
given by [22] 


Ck = 



(52) 


Both dehnitions (50) and (51) are valid for a free held system in equilibrium, and have proven to 
be very useful in interacting theories out of equilibrium as well [21, 30, 68, 69]. From the studies in 
1-|-1 and 2-|-l dimensions [48, 49], we expect the system to exhibit a quasiparticle structure before 
reaching thermal equilibrium. The dehnitions (50) and (51) can be used to monitor the evolution 
of the system towards such a quasiparticle-like state, and eventually to equilibrium. 

Once the system is close to equilibrium, we can read from (51) the effective quasiparticle mass 
mes{t) by comparing it to the dispersion relation 

ujlit) = c^{t) [mlffit) + k^) , (53) 


where the factor c{t) is a measure of an effective speed of light or an inverse refractive index. 
A temperature Teff(t) and chemical potential Hef[{t) can be determined by htting the occupation 
number (50) to a Bose-Einstein distribution 


np{t) 


1 


(54) 


In 


1 

1 H- 

Up 


1 



using 


heS 

TeS 


(55) 
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We also keep track of the ‘memory kernels’ in the equations of motion (29-32), i.e. the self¬ 
energies t'), t') and t'), which can be compared with perturbative estimates. Limits 

on computer resources (memory and CPU time) requires us to cut the memory kernels and thus 
keep only some finite range backwards in time (i.e. 0 for \t — t'\ > tent)- The size of the 

self-energies helps us determine whether the cut was late enough for the discarded memory integrals 
to be negligible. A way to judge whether the discarded memory was indeed unimportant for the 
dynamics is to verify that the total energy density (35) is conserved. Monitoring the evolution of 
the energy density (35), one finds that at very late times, the effect of the memory cut shows up 
as a very slow drift in the energy. In all the runs presented here, the energy is conserved to within 
2%. For smaller lattices we checked that later memory cuts make the drift smaller. We found no 
such drift if the whole kernel was kept. 


D. Symmetric phase: equilibration 


We hrst consider the evolution of the system in the symmetric phase. The simulations are 
performed on a lattice with = 16^ sites, lattice spacing am = 0.7, time-step dt = 0.1 and 
coupling^ A = 6. The memory kernel is cut off at mfeut = 28 unless otherwise specified. In the 
following, all quantities shall be expressed in renormalized mass units, i.e. units of m. 

For the initial conditions of the propagators we shall take: 


• Thermal, with Tin/m = 1.36,1.43,1.93,2.86. 

• Top-hat 1 (Tl), with = 2.04, = 6.12 and t] = 2. 

• Top-hat 2 (T2), with k^jg/m^ = 0, k^g^/m^ = 5.71 and rj = 1.85. 

• Top-hat 3 (T3), with k^jg/m^ = 6.12, k^g^/m^ = 8.16 and rj = 1.6. 

The three top-hat initial conditions have the same initial energy. 

In a quasi-particle picture, we can introduce the total particle number density 

Wot f f for infinite volume in the continuum, 

^tot = 1 1 11- 

^ 17^ Wk on the lattice. 


(56) 


The three top-hat initial conditions T1-T3 do not have the same total number of particles, although 
Tl and T2 are fairly close to each other. More about this below. 

For the initial mean field we take (/>(0) = 0. In this case, the Hartree and the two-loop ap¬ 
proximations are identical. In the following we consider the evolution in both the Hartree and the 
“basketball” approximations. 

In Fig. 2, we show the evolution of versus uik, starting from the Tl initial condition. The 
Hartree approximation (black) is compared with the “basketball” approximation (green/grey). In 
the former case, there is no equilibration. For the “basketball” case, we observe that the energy in 
the excited modes is distributed via scattering. As we shall see this leads eventually to a thermal 
distribution. 

In Fig. 3 we follow the evolution of the dispersion relation with the same initial condition, again 
comparing Hartree to “basketball”. Notice the oscillating pattern early on in both cases. In the 
“basketball” approximation the modes eventually relax to a perfect straight line. It turns out that 


Recall that we neglect the difference between A and Aq. 
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Figure 2: Evolution in time of the occupation numbers rik vs. Wk, for the T1 initial condition. We display 
the results of the “basketball” (black dots) and the Hartree approximation (green/grey dots). 


at the couplings and energies used here, the coefficient cP{t) is equal to 1 up to well within one 
percent. We will therefore assume it to be 1 in the following. Although we do not show it here, 
we found that larger coupling and large energy density results in a faster evolution towards this 
quasi-particle state. 

Judging by eye. Figs. 2 and 3 suggest that already at times mt = 56 to 84, the system behaves 
as aproximately thermal, for T1 initial conditions. Still, this is presumably much later than a 
pre-thermalization time based on the the equation of state, as studied in [4]. However, it does not 
mean that the memory of the initial conditions in the particle distribution is already lost by times 
mt > 84. 

It was remarked already in the studies in 1-|-1 [44] and 2-|-l dimensions [49] that the final state 
depends only on the energy density (at a given coupling). Indeed, it was found that the limit 
distribution function nk corresponds to a Bose-Einstein, characterized by just one parameter, the 
temperature. Fig. 4 shows the evolution of individual modes when starting from the Tl, T2 and 
T3 initial conditions, which have the same energy density. For Tl and T2 we see that the modes 
approach a common final value. It seems reasonable to call this stage kinetic equilibration, as the 
kinetic energy is redistributed over the modes to reach a Bose-Einstein distribution. However, as 
we will see below, the total number of particles is not adjusting as fast and it still remembers the 
initial state by the time kinetic equilibration is completed. 

As mentioned before, the initial condition T3 has not only a different initial spectrum, but 
also a different total number of particles. It also reaches kinetic equilibration, but with a different 
kinetically equilibrated state. 

At intermediate times {mt ~ 1000) kinetic equilibration has taken place, and we can compare 
the distribution functions and dispersion relations, Fig. 5. Tl and T2 have equilibrated to al¬ 
most identical Bose-Einstein distribution functions, parametrized by an effective mass, an effective 
temperature and an effective chemical potential. T3 has reached a different Bose-Einstein with 
a different temperature and chemical potential, and a slightly different effective mass. We have 
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Figure 3: Time evolution of the dispersion relation vs. k^) in the symmetric phase, starting from the 
T1 initial conditions, in the “basketball” ( black dots) and Hartree approximations (green/grey dots). 



Figure 4: Evolution of individual modes for a T1 


T2 and T3 initial condition with same energy density. 


included a number of thermal initial conditions for comparison. By construction, these have no 
initial chemical potential and remain so to a very good approximation. 

Whereas kinetic equilibration can be the result of simple 2 2 scattering, chemical equilibra¬ 
tion, which changes the total particle number, happens through 1 3, 2 4 and higher order 

processes. These are included due to the resummations performed by the <l>-derivable approxi¬ 
mation into the “sunset” self-energy diagram. Approaches that only take into account on-shell 
scattering, such as the Boltzmann equation with only binary collisions 2 2, cannot account for 

chemical equilibration. What we see is that kinetic equilibration including memory loss happens 
on a time scale of about 500 — 1000/m, whereas chemical equilibration is a much slower process. 
Effectively, there is a chemical potential in the initial stages, causing initial conditions with different 
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A^tot to relax to different intermediate kinetically equilibrated states. 



Figure 5: Left: occupation numbers ln(l + 1/uk) at intermediate times mt = 1000 for initial Tl, T2 and T3 
and a set of thermal initial states. Right: The dispersion relation for the same cases. The intercepts give 
the effective masses squared. 


We illustrate this point in Fig. 6 , left-hand plot, where we show the evolution iVtot for the 
initial conditions T1-T3. For comparison, we also include the Wot of Bose-Einstein distributions 
at various temperatures. In the right-hand plot we follow the evolution of the effective mass, 
temperature and chemical potential for the Tl case. The time evolution can be well reproduced 
by exponential fits of the form a* -|- bi exp {—'jit) (the dashed lines in the plot), suggesting an 
asymptotic temperature of around T/m = 1.36. Also, within a factor of two, fits to the three 
quantities all suggest an equilibration time of around 7 ^”^ ~ 10^/m. Chemical equilibration is a 
full order of magnitude slower than kinetic equilibration in this system. Comparing with the study 
in [49], it appears that chemical equilibration is much slower in 3-1-1 than in 2-|-l dimensions. The 
fit to the chemical potential is not as good as to the effective temperature or mass, and it also 
predicts a non-zero asymptotic chemical potential (/r/m = 0.7). This is consistent with our above- 
mentioned interpretation that the system is in a prethermalized stage [4] for which an exponential 
extrapolation of the evolution of T and ^ does not necessarily yield the actual asymptotic values. 

As an aside we compare the observed mass with an estimate that results from the Hartree 
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Figure 6: Left: Total particle number density ntot versus time for the initial conditions Tl, T2 and T3. 
The dotted lines are ritot for Bose-Einstein distributions at different temperatures. Right: Evolution of 
the effective mass (squared), temperature and chemical potential for late times starting from the Tl initial 
condition, with exponential fits (dashed lines). The Hartree estimate for the mass, (31), and the solution of 
the gap equation (60), are also shown. 
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approximation. At a given time t, the finite gap equation for the Hartree effective mass reads 

A C 

=m^ + T}{t) - 5rr? = rr? + - J t) - Sm^. (57) 

The mass counterterm 5m? is given by the vacuum part of the “leaf” self-energy diagram, as 
described previously. For the correlator F{t,t) in Eq. (57) we take the same form as a free quasi¬ 
particle gas in equilibrium, i.e. 




1 

UJk{t) 


nUt) + - 


(58) 


Here nk(t) is a Bose-Einstein distribution function with the temperature and chemical potential 
obtained from the simulations at time t, and uj]s_{t) is here dehned in terms of the effective mass as 

The result for the Hartree effective mass Mh is then determined by the 
self-consistent gap equation 


fi) — m — 


2 X f d?k n]i{t;T,fi) 2 


(27r)3 a;k(t) 


= m + 


XT^ 


dx 




(59) 


To compare with the numerical result we use the lattice analog of the gap equation (59), i.e. 


n) — m — 


A 1 


2 (Na) 


E 


nk(t;r, n) 


(60) 


For the case of the evolution starting from the T1 initial condition, the lattice Hartree mass is 
shown in Fig. 6. We see that it is slightly higher than the effective mass obtained in the simulation 
with the “basketball” approximation. At least in this case, the contribution from the “sunset” 
diagram to the mass appears to be small relative to the Hartree case. 


E. Symmetric phase: Damping and the spectral function 

1. Mean field damping 


We now consider a situation already in (or close to) thermal equilibrium, where the mean field 
is slightly displaced from its equilibrium value (/> = 0. This allows us to study the response of the 
system to small perturbations. In this case, the mean field evolution can be studied by linearizing 
the equation of motion (32) around the equilibrium value. For homogeneous fields this leads to 

0(t) + m2 (T, t)(t>{t) + [ dt' S(; (t, t') (fit') = 0. (61) 

do 

with SQ(t, t') the zero momentum mode of the “sunset” self-energy and M (T, t) given by (31). Close 
enough to equilibrium we may assume time translation invariance, such that t') depends only 
on t — F and M{T,t) is constant. The equation (61) can then be solved by a Laplace transform in 
the time coordinate [70]. Taking as initial conditions for the mean field i;A(0) = (pi and 0(0) = 0, 
the solution to the linearized equation of motion (61) can be written as [27, 70] 


(p{t) 



ct;ImSQ(a;) cos(a;t) 
M2 - ReS^(u;) 


+ ImS^(u;)2 


(62) 
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where ReS^ and ImS^ correspond respectively to the real and imaginary part of the retarded 
self-energy, given by y) = 0(xo — yo)f^^{x, y). For weak coupling there is a narrow resonance 

at (u = Meff, with -|- ReSg ('^)- To a good approximation, one finds that for short times 

the evolution is given by [71] 


with 


(/>(t) (l)iZ e cos (Mefft — a ), 


Z = 


a = 


1 - 


aReS^(Meff) 




-1 


dlmt§{M,s) 




7 = Z 


ImSg(Meff) 

MeS 


(63) 


(64) 

(65) 

( 66 ) 


The parameter 7 corresponds to the on-shell damping rate. For weak enough couplings one can 
approximate Z ~ 1 and Mgfr = M for the calculation of 7 . From ( 66 ) we see that the damping rate 
is determined by the imaginary part of which corresponds to the “sunset” self-energy diagram. 
In the context of perturbation theory, Im can be calculated analytically and the damping rate 
is found to be [72] 


7(M) 


1287r3M 


Li 2 (e-^/^) 


(67) 


where \j\ 2 {z) is the second polylogarithmic function, defined by Spence’s integral 


U2{z)=- r 

Jo 


dw 


In (1 — w) 


w 


( 68 ) 


For temperatures T ^ m, the damping rate follows from the expression of the high-temperature 
screening mass [73] 


+ 


Ar2 A 


__Mr + 0 AMln— . 


24 Stt 


m 2 A 


T2 J 


(69) 


In the limit of very weak coupling and high temperature one obtains for the damping rate the 
compact result^ [72, 74] 


7 


T>m, A<1 


A 2 r 2 

7687rM ■ 


(70) 


In the numerical simulations the mean field is initially displaced to the value (pi/m = 0.142. 
For such a small perturbation, the mean field is expected to perform a damped oscillation of the 
type (63). We fit the evolution of the mean field with Eq. (63), which allows us to extract the 
effective mass Meff and the damping rate 7 . These will depend on the strength of the coupling 
and the temperature. The behavior of the mean field is studied in a thermal bath at temperatures 
T'm/'m = 0,1.43,2.86 for both the two-loop and “basketball” approximations. For the “basketball” 


® This approximation for the damping rate is often used in the literature. For the values of the coupling A = 1,6 
used in the numerical analysis presented in this paper, however, the approximation (70) to (67) is not valid, even 
for high temperatures. 
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Figure 7: Time evolution (pit) of the initially displaced mean field (pi = 0.142 for the cases Ty^jm = 
0,1.43,2.86, in the “basketball” approximation. 


1.25 

1.2 



' 1 ' 1 

/. 

/■ 

k ' 

0.02 

- ^^^^^ 

1 ^ r 


• • Haitree (Continuum) 




— • Perturbative (Continuum) 

/ 


— • Haitree (Lattice) 

/ 




■ ■ Perturbative (Lattice) 

/ A 


A Zero mode mass, Basketball 

/ 



- 

■ Mean field damping, Two-loop 

/ 


■ Mean field mass, Two-loop 

/ ♦ 




A Zero mode damping, Basketball 



• Mean field mass. Basketball 

/ 



• 

• Mean field damping, Basketball 

/ I 



/ 

/ • 


0.01 

- 


/ ■ .’ 

/ 

- 


/ 

/ 

/ * 


S 




/ 

/ 


u 

✓ 



/ 

< 







/* 







y 



o' 

L 


- 








ir 

^^_ 

^^ _L 



1_1_1_1_i 

1 1 1 1 


0.5 


1 1.5 

T/m 


0.5 


1 1.5 

T/m 


Figure 8: Left: Effective masses from the analysis of the evolution of the mean field (in the two-loop and 
“basketball” aproximations) and spectral function zero-mode (“basketball” case). Hartree estimates, both 
in the continuum (dotted line) and on the lattice (dashed line), are also shown. Right: Damping rates. 


case, the mean field evolution is shown in Fig. 7. As mentioned earlier, the damping of the mean 
field is present in both the two-loop and “basketball” approximations. The difference between the 
two cases lies in the fact that the correlators F and p evolve quite differently. For the two-loop 
case, the equations of motion for F and p contain almost no damping, since the only potential 
contribution to damping is in the “eye” diagram, which is proportional to (p"^, and thus tiny for 
(/> ~ 0. For the “basketball” case, however, the equations of motion for F and p contain damping 
through the “sunset” diagram. The differences in the evolution of the correlators for the two-loop 
and the “basketball” approximations enters as a higher-order effect in the evolution of the mean 
field. In particular, this may lead to different effective masses and mean field damping rates. 
We show these differences for various temperatures in Fig. 8, where the results for the two-loop 
(squares) and “basketball” (large dots) are plotted. 

For comparison we evaluated the perturbative result (67), using the Hartree mass (59) for M. 
To see the finite-volume and discretization effects we also did the analytical computation on a 
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spatial lattice. The Hartree mass is in this case given by (60). In finite volume, the discreteness of 
the momenta leads to complications in the calculation of the damping rate, which we dealt with 
along the lines presented in [75]^. The perturbative results for the mass and damping rate are also 
presented in Fig. 8. As we can see from the mass plot (left), the correction to the mass coming from 
the “basketball” approximation is small relative to the Hartree case. In the damping plot (right) 
we observe that the damping rate obtained from the numerical analysis of both approximations is 
substantially larger (about 20-40%) than the perturbative result (on the lattice). The continuum 
and the lattice perturbative results begin to differ around T/m > 1 due to cut-off effects. 


2. Propagator damping and spectral function 


Damping in the propagator can be elegantly phrased in terms of the spectral function pk^(t,t'). 
In a situation close to thermal equilibrium we expect it to be time translation invariant and in a 
narrow-width approximation be given by 

Pk(i, ^0 = —sin [wk(t - t')] > (71) 

To study the approach to equilibrium of the spectral funcion, it is useful to perform a Wigner 
transformation in terms of the mean time T = {t +1')/2 and relative time t = t — t'. This can be 
written as 

f2r 

p]^{uj,T) = 2i dTsm{u}T)pk{T + t/2,T - t/2). (72) 

Jo 

Since we are solving the equations of motion in a finite time and keep information only as far back 
as the memory kernel, we have a cut-off in the integral of (72) as 


cAut 

Pk{u},T) ^ pi,{uj,T;tcnt) = “Ji dTsm{uJT)pk{T + t/2,T - t/2). (73) 

Jo 

If the system is sufficiently close to equilibrium, time translation invariance should be a good 
approximation, p\i{t,t') ^ pk{t — F), and 


r Aut 

/)k(w,7';tcut) ~ 2z / dr sin (ujt) pk{T,T - t). (74) 

Jo 

With this approximation the integrand in (74) runs from p\^{T,T) to p]i(T,T — tent), which is 
convenient for numerical purposes. In the following we shall make use of this approximate Wigner 
transform. 

The approximation (74) is valid provided T is large enough so that the system is close to thermal 
equilibrium. In that case Pk{(^\i,'T) should be well approximated by a Breit-Wigner form 




4a;(7k/2) 

(u;2_^2)2+^2(^^/2)2- 


(75) 


The evolution of the spectral function pk.{t,t') starting from a thermal background at T/m = 2.86 
and A = 6 is shown for two different kernel lengths, mfeut = 28 (green/grey) and mtcut = 84 


^ For example, equations like (101) in the appendix do not make sense anymore. We evaluated the frequency 
integral in the solution for the linearized equation of motion (62) for fit), using a finite “ie” in the retarded sunset 
self-energy on the lattice. 
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Figure 9: Left: The evolution of the spectral function po{t, t') starting from a thermal initial condition with 
T/m = 2.86 and with coupling A = 6. The result with a short and a long kernel is plotted. Overlaid a fit to 
the form (71). Right: The approximate Wigner transform po{uj) at time mT = 200. A Breit-Wigner fit to 
the Wigner transform with kernel 84/m is included (dashed line). 


(black), in Fig. 9 (left). Overlaid, although barely discernible, a fit of the form (71). As can be 
seen, the fit is excellent and it is the same fit for the two kernel lengths. In Fig. 9 (right) we show 
the result of the approximate Wigner transforms for the spectral function of the spatial zero mode 
at mT = 200 and with two different kernel lengths. For the long kernel case, one can nicely fit a 
Breit-Wigner form, which gives the same damping rates and masses as in the fit of the left-hand 
plot. For the short kernel case, the Breit-Wigner fit is not so accurate. In the following analysis, we 
extract the damping rates and masses from fits directly to the time-representation of the spectral 
function with kernel length mtcut = 28. 

Fig. 8 shows the dependence on temperature of the effective mass (left) and damping rate (right) 
from fits of the form (71) for the spectral function zero-mode together with the fits (63) 

for the mean field discussed earlier. For pQ{t,t'), we show only the results for the case of the 
“basketball” approximation, since, for the values of the mean field considered here, it is practically 
zero in the two-loop approximation. The spectral-function zero-mode mass and damping rate 
(plotted with triangles) closely follow the values for the mean field. 


F. Broken phase: equilibration 

A similar analysis can be carried out in the broken phase, where there is a non-zero mean field 
present. In this case we can use both the two-loop and “basketball” approximations to study the 
damping of the correlators. From the point of view of perturbation theory, there is no damping 
in the two-loop approximation, for which only the perturbative “leaf” and the “eye” diagrams 
contribute to the self-energy, and their imaginary parts vanish on-shell (see also appendix B). Our 
task will be to study the damping in the two-loop approximation from the <I>-derived equations of 
motion. These formally take into account the contributions from all orders in perturbation theory 
that result from any iterated insertion of the “leaf” and “eye” diagrams into the self-energy. These 
diagrams contain off-shell scattering effects that can, in principle, lead to a total non-zero on- 
shell imaginary part for the self-energy, thus providing damping. For the case of the “basketball” 
approximation, the “sunset” diagram enters in the self-energy (13), which contributes to damping 
even in perturbation theory [72, 74]. The solution of the <I>-derived equations of motion leads to 
additional contributions from higher orders compared to perturbation theory, and thus one expects 
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Figure 10: Evolution of the occupation numbers (rik vs. LUk_/m), starting from the T1 initial condition, 
in the broken phase (A = 1) and in the two-loop (green/grey) and “basketball” (black) approximations, 
respectively. 


to find a larger damping and faster equilibration. 

Approximations based on truncations of the loop expansion of the 2PI effective action suffer 
from instabilities which make it impossible to treat very large couplings and very large energy 
densities or particle numbers. In this sense, the <I>-derived equations of motion can be thought of as 
resummed perturbative, useful in the domain of weak coupling and small fields. In the symmetric- 
phase simulations described previously, A = 6 is in the upper end of what stays stable in our 
experience, whereas we can use temperatures (or energy densities corresponding to temperatures) 
up to T/m = 6 or even higher. In the broken phase, the instabilities turn out to be even more 
constraining. In particular, we will need to use a smaller coupling (A = 1) and temperatures 
below T/m = 2. As we have seen, the latter is not much of a restriction since it still covers the 
region where cut-off effects are small. However, it implies that equilibration and damping takes 
much longer (the damping times scale roughly as A^T^ for the sunset diagram). In particular 
we need to use a longer time kernel (we use mtcut = 84) and we will not be able to track the 
evolution far enough to see chemical equilibration. We shall content ourselves with establishing 
kinetic equilibration and studying the damping of the mean field and the modes of the spectral 
function. We have no doubt that chemical equilibration will take place as well. In particular, we 
will see that total particle number is not conserved. 

The mean field is taken initially to be at the tree-level value = 0) = utree- This is not the 
self-consistent finite temperature solution of the truncated equations of motion, but a bit displaced 
from it. Due to this initial displacement, the mean field will oscillate and damp to its equilibrium 
value. For the propagators we will use thermal initial conditions, as well as top-hat 1 (Tl). Notice 
that the input mass is the broken phase zero temperature mass \/2m rather than the symmetric 
phase one. All results are still in units of m. 

The evolution of the occupation numbers and the dispersion relation for both the two-loop and 
“basketball” approximations are shown in Fig. 10 and Fig. 11. Both cases show that (kinetic) 
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Figure 11: Evolution of the dispersion relation vs. kF'/m^), with the T1 initial condition, in the 

two-loop (green/grey) and basketball (black) approximations, respectively. 



Figure 12: Left: The total particle number density ntot(t) for the two truncations, normalized to the initial 
value. Right: Particle distributions at the latest time mt = 1000, for T1 and thermal initial conditions. 


equilibration is taking place. In the “basketball” case, equilibration is slightly faster. Interestingly, 
the off-shell scattering effects taken into account by the 2PI effective action with only the eye 
diagram lead to an equilibration almost as fast as in the “basketball” case. Chemical equilibration 
happens on much longer time scales, and although we found that the total particle number does 
change in time (Fig. 12, left), the reach of our simulations was insufficient to estimate the asymptotic 
temperature. At our latest time of mt = 1000, the distribution is consistent with a Bose-Einstein 
with T/m = 1.24 and /i/m = 1.12 (Fig. 12, right). 
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Figure 13: Left: Evolution of the mean field in the broken phase in a thermal background at T/m = 2.14, 
in the two-loop and “basketball” approximations. Right: Effective masses, vacuum expectation values and 
damping rates from fits of the form (63) for both the two-loop (empty symbols) and “basketball” (filled 
symbols) approximations. The lines represent the Hartree approximation. 


G. Broken phase: Mean field damping and the spectral function 

For weak coupling we expect the position of the equilibrium mean field expectation value v to 
be close to the initial value (/)(0) = Utree- In the case of thermal initial conditions, the initial mean 
field displacement corresponds to a small perturbation. As in the symmetric phase case, one can 
study the evolution of the mean field by linearizing the equation of motion around the equilibrium 
value. We write (/)(t) = v + a{t), where a{t) is the deviation. The linearized equation of motion for 
a is then 

a{t) + M‘^{T,t)a{t) + f dt'= 0. (76) 

Jo 

The vacuum expectation value v is the solution of 

\ ft 

M^{T,t)v - -v^ + / = 0. (77) 

3 Jo 

For weak coupling 


(78) 

In Eqns.(76-78), M{T,t) corresponds to the finite temperature Hartree effective mass in the broken 
phase, given by 

= -rr? + ^ J 

The analysis of the evolution of a proceeds as in the case of the symmetric phase. For weak enough 
coupling the mean field damping rate is approximately given by the perturbative estimate (67). 
Fig. 13 (left) shows the mean field evolution at T/m = 2.14, in the two-loop and “basketball” 
approximations, respectively. In both cases, the mean field performs a damped oscillation, from 
which we can extract an effective mass and mean field expectation value. As we can see, the 
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Figure 14: Wigner transforms of the spectral function for the four lowest modes, starting from a T1 initial 
condition, at two times (left plot, mT = 252 and right plot, mT = 1000) for the “basketball” (full lines) 
and two-loop (dotted lines) approximations. The weak T-dependence of the former shows that the system 
is close to equilibrium. Indeed, Breit-Wigner fits work nicely (dashed lines). The strong T-dependence in 
the two-loop case suggests that the system is not yet sufficiently close to equilibrium. 


damping does not follow a simple exponential form, and curiously, the two-loop data appear to 
indicate faster damping than the basketball data. Still, we shall use an exponential fit as a rough 
estimate of the damping rate. The temperature dependence of these frequencies and damping 
rates, in addition to the field expectation values, is shown in Fig. 13 (right) for the two-loop (open 
symbols) and “basketball” (filled symbols) approximations. The masses and field expectation 
values are indistinguishable for the two truncations. They are slightly off the respective Hartree 
estimates (77) and (79) (full lines), indicating as in the previous section that the contribution of 
the “sunset” diagram in the “basketball” approximation is small relatively to the Hartree case. 
The damping rates in the two approximations are consistent with each other. 

Similarly, we observed damping in the evolution of the spectral function pk(^) t') in both approx¬ 
imations. This damping is small and not well approximated by an exponential form. We performed 
the Wigner transformation as specified in Eq. (74), see Fig. 14 to the data of the late time evolu¬ 
tion starting from the T1 initial condition. The value of the cut-off tcut = 84 produces some 
noise, but for the “basketball” approximation (full lines) there is clearly a well determined peak 
with a finite width at all times, which can be fit with a Breit-Wigner form (dashed lines). In the 
basketball approximation, the form of the Wigner transform does not change much in time, which 
indicates that the system is relatively close to equilibrium. However, the results are different for 
the approximate Wigner transforms in the two-loop approximation (dotted lines). In this case 
the transformed spectral function changes significantly in time from thin peaks early on to less 
clear maxima at later times, still localized around the peaks, see Fig 14 (right). We do not fully 
understand the reason for this discrepancy; it could be, that the cut-off time fcut = 84 in our 
implementation of the Wigner transform is too short in case of the two-loop approximation. On 
the other hand, the large T-dependence of the spectral function indicates that the system may 
not be sufficiently close to equilibrium in the two-loop approximation, in which case the use of 
the approximate Wigner transform (74) is questionable. The approximation itself appears to work 
reasonably well, as shown in Figs. 12 and 13. 
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V. CONCLUSIONS 

We have studied equilibration in theory in 3+1 dimensions for a variety of initial conditions, 
both in the symmetric and broken phase. Two different <l>-derivable approximations including scat¬ 
tering effects have been used: two-loop and “basketball”, the latter corresponding to a truncation 
of the 2PI effective action at O(A^). In the symmetric phase the two-loop and the “basketball” 
approximations differ in that the first includes damping into the evolution of the mean field only, 
whereas in the second it is also present in the equation of motion for the 2-point functions. In the 
broken phase both approximations include scattering effects into the 2-point functions and thus 
can lead to equilibration. 

From the numerical study of the evolution of the occupation numbers we were able to establish 
that in the symmetric phase both kinetic and chemical equilibration is taking place, the latter at 
a substantially slower rate. 

By analyzing various initial conditions we found that, after kinetic equilibration, the occupation 
numbers at intermediate times are given by a Bose-Einstein distribution with an effective chemical 
potential. This is similar to what was found in previous studies in 2+1 dimensions [49]. Given the 
same total energy, we found the intermediate effective chemical potential to be generally non-zero, 
with its size related to the initial total particle numbers, as may be expected. Eventually, one would 
expect the limit distribution to depend only on the energy density of the system, and close to a 
Bose-Einstein with zero chemical potential [44, 49]. Comparing to the studies in 2+1 dimensions 
[49], our numerical analysis indicates that the subsequent chemical equilibration is much slower in 
3+1 dimensions. 

We were also able to extract effective masses and damping rates from the analysis of the evolu¬ 
tion of the mean field and the spectral function. The contributions to the mass from the two-loop 
and the “basketball” approximations seem to be small comparing to the Hartree case. In the sym¬ 
metric phase, the results for the damping rate are about a 20-40% higher than the perturbative 
estimates. This indicates that the scattering effects associated with the resummations encoded in 
the <h-derivable approximation are substantial. Finally, we checked that the damping rate obtained 
from the mean field coincides with the one from the spectral function zero-mode. 

In the broken phase we found that both the two-loop and the “basketball” approximation lead 
to equilibration. Surprisingly, the equilibration seems to be just a bit slower in the two-loop case. 
This is particularly remarkable, since, in perturbation theory, the two-loop approximation does 
not have on-shell damping. Indeed, in that case only the perturbative “leaf” and “eye” self-energy 
diagrams contribute, and their imaginary parts vanish on-shell (see appendix B). The fact that the 
two-loop approximation in theory equilibrates so fast might be relevant to pure gauge theories, 
where the lowest order ^-derivable approximation (at 0{g^)) considers the same diagrams. 

On a practical note, we found that the loop expansion suffers from restrictions reminiscent 
of perturbative expansions, in that large couplings and/or large field occupation numbers trigger 
instabilities when solving the equations of motion. This has to our knowledge not been reported 
for simulations in I+I and 2+1 dimensions, although we have found them in those cases as well. 
In addition to the instabilities, issues such as CPU time and computer memory necessary for 
dealing with the memory integrals are significant restrictions, especially when studying late time 
thermalization. Expansions in \/N with N the number of fields have been shown to be more stable 
and able to deal with non-perturbatively large occupation numbers at large coupling [31, 52, 54, 55]. 
In such cases care should be taken to ensure that the expansion is controlled by using a sufficiently 
large value of [41]. 
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APPENDIX A: ENERGY-MOMENTUM TENSOR IN 4>-DERIVABLE 

APPROXIMATIONS 


The energy-momentum tensor for a given truncation of the 2PI effective action can be deter¬ 
mined through Noether’s procedure, i.e. by identifying the current term resulting from the space- 
time dependent translations ^ + e^{x). A convenient way to find the Noether current is 

to pertorm an infinitesimal translation that vanishes on the space-time boundary. The translation 
x^ ^ x^ + €^{x) can be viewed as a transformation of the relevant variables, which in the case of 
the 2PI effective action are the mean field 4>{x) and the connected 2-point function G{x,y). This 
transformation is given by 

(/)(x) —> = (j){x + e(x)) = fi{x) + €^{x)df^fi{x), (80) 

G(x, y) —> G'{x, y) = G{x + e{x),y + e{y)) = G{x, y) + €^(x)d^G(x, y) + e'^(y)5f G(x, y), (81) 

where the variables that the partial derivatives act on are indicated with a superscript. Under 
these transformations the variation of the 2PI effective action r[i^, G] can be formally written as 

6T[fi,G] = [ T^"(x)a^e,(x), (82) 

J X 

where f^ = f d^x. The quantity T^‘' defines a conserved Noether current, which is identified as 
the energy-momentum tensor. To see that it is indeed conserved, notice that when the <I>-derived 
equations of motion for and G are satisfied. This applies, in particular, to the transformations 
(80) and (81), and hence (ir[()), G] = 0. Taking (82) and making a partial integration one obtains 


5r[fi, G]= - f e,(x) a^r^"(x) = 0. (83) 

X 

Since ei^{x) can be taken arbitrary, the energy-momentum tensor is conserved, i.e. d^T^'^{x) = 0. 
Below we give the explicit form of the energy-momentum tensor for any d>-derivable approxima¬ 
tion by applying the transformations (80-81) and using the definition (82). We study independently 
the contributions coming from the four terms in the action (2): 


i) 


The first term in Eq. (2) is given by ^[c))] and leads to the usual form of the energy-momentum 
tensor for the mean field (j), namely 


rf^(x) = a^<)>(x)aXx) - 


^dx(l){x)d^cl){x) - ^m^cj){xf - 


(84) 


ii) The second term appearing in Eq. (2) does not contribute to the energy-momentum tensor. 
Indeed, applying the transformations (81) leads to 


(5[TrlnG] = TVG-^<5G 

G~^{y,x) \e^{x)dxG{x,y) + {x ^ y) 


X jy 


= 2 


’(x)9f(l*-"^^(x — x) = 0. (85) 


Hence, the corresponding contribution T^^(x) to the energy-momentum tensor vanishes. 
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iii) To obtain the contribution from the third term, given by (i/2) Tr [Gq ^ — G G, one pro¬ 
ceeds similarly to what was done in i). Under the transformation (81), it becomes 

5 ItV [Go'- G-i] G = 

^ J J -m^ - {x - y) e^{x)dlG{x, y) + {x ^ y) 

-y/ J (t>{x)e^{x)dl(t){x)5^'^\x-y)G{x,y). (86) 

After some partial integrations and making use of the identity 

[ (x - y)G{x, y)l = [ (5^^^ (x - y) [d^G{x, y) + d^G{y, x)] , (87) 


one can write the above in the form given by Eq. (82), which allows to extract the contribution 
of this term to the energy-momentum tensor. One finds 

T,^^{x) = I 5 ^^Hx - y) + y^^A<)>(x)2j G(x,y). (88) 


iv) For the fourth term in Eq. 2, which is given by the functional $[(/>, G], the transformations 
(80-81) give 




6G{x,y) 


ef,{x)d^G{y,x) + ef,{y)d^G{y,x)] . (89) 


What we want is to write this in a form similar to (82) such that its contribution to the 
energy-momentum tensor can be extracted. To do this, notice that the functional G] 
is a scalar quantity that does not contain derivative terms. This means that, under the 
space-time translation x^ ^ x^ -|- e^(x), the terms in $ only change by the appearance 
of the Jacobian of the transformation at every loop integration. This Jacobian can be 
accommodated by a simultaneous change in a scale factor ()(x) introduced at every integration 
vertex as A ^ C(3;)A [46, 76]. Thus the simultaneous variation 

(/(x) ^ /)(x-be(x)), G(x, y) G{x+e{x),y+e{y)), ({x) = 1 ^ C(x) = det {6^ + d^e^{x)) 

(90) 

leaves the functional d* invariant. For infinitesimal transformations, this invariance implies 

I + ^Ay)BSO{y.x)]+l = o. 

(91) 

One can then use the identity (91) to write the variation JJ>[i/), G] in a form similar to (82) 




y.jc(x)ic=i ^ ^" 

In this manner, the contribution of the functional <h to the energy momentum tensor can be 
written as 


Ttix) = 


SC(x) lc=i‘ 
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The total energy-momentum tensor is obtained by adding up all the contributions from i)-iv), 
i.e. T^^(x) = T|‘^(a;) -|- T!^'^{x) + -|- T!^'^{x). The result can be compactly written as 


r^^(x) = 


(0(x)0(2/) + G(x, y)) 


x=y 

2 / 




+ g^'^^Xct^ixf + ig^>(x)^G(x, x) - 


C=i 


(94) 


APPENDIX B: PERTURBATIVE DAMPING FROM THE ’’EYE” DIAGRAM 


For completeness, we include here the calculation of the imaginary part of the perturbative 
“eye” diagram in equilibrium in the real-time formalism, using the Schwinger-Keldysh contour. 
More details can be found in e.g. [77]. We introduce the labels -|- or — to indicate whether the time 
variables of any quantity live respectively on the C~^ or C~ branch of the contour. In terms of the 
various contour components, the retarded self-energy is given by S^(x,y) = S''''''(x, y)-l-S'' (x,y). 

For the case of the eye-diagram, in momentum space one has 

Eg(p) = S++(p) + S+-(p) = ^ [G++(k)G++(k -p)- G+-(k)G-+(k - p)] , (95) 

where v{T) is the mean field equilibrium expectation value at temperature T. We shall use 
and to denote the 4- and 3-dimensional momentum integrations f d‘^k/{2'K)^ and f d^k/{27r)^ 
respectively. 

It is convenient to use the Keldysh basis [10, 13], where the various components of G are given 
in terms of the symmetric, retarded and advanced correlators F, G^ and G^ respectively. Their 
perturbative expressions are 


F{k) = 2TT5{k^ — w?) n{kQ) + X 


G^{k) = 

a^(k) = 


1 


(/cq + FY — — m? ’ 

1 

{ko — ieY — ’ 


with e = 0"*“ and n(A:o) the Bose-Einstein distribution at temperature T and energy ko. 
In the Keldysh basis the retarded self-energy (95) becomes 


(96) 

(97) 

(98) 


S«(p) 



{GR{k)GR{k -p) + GA{k)GA{k - p))-iF{k)GR{k-p)-iGA{k)F{k-p) 


(99) 

The first two terms in the RHS have poles at only one side of the complex plane. In the integration 
over /co one can always choose to close the contour at the other side, thus these two terms vanish. 
The last two contributions can be seen to be equal to each other by the change of variable k 
(p-k). After performing the /cq integration with the help of the 5-functions in F we obtain 


S§.(a;, p) = \^v{Tf + 


1 

2 

+ 


1 


(cu - Wk)^ - - iesgn{uj - Uk) 

1 

(w -h UkY - ‘^p-k “ sgn((u -F (Uk) 


( 100 ) 
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The imaginary part of the self-energy is obtained by using l/{x + ie) = P{l/x) — iTT6{x) and 
decomposing the resulting delta functions. For a; > 0, and after convenient changes of variable, we 
obtain 


ImS:^(a;,p) 


A 2 n(T )2 

2 


/ " 

_/k k 


{[2nk -h 1] S{uj - cuk - a;p„k) + 2n^6{uj Wk - Wp-k)} • 


( 101 ) 

The first contribution to the integral corresponds to the decay of an off-shell excitation into two 
on-shell excitations. The second one corresponds to Landau damping via scattering of the off-shell 
excitation with on-shell particles from the heat bath (occuring only at T 7 ^ 0 ). 

One can make use of the delta functions present in (101) to solve the angular part of the integral 
over the internal momentum k. Indeed, using the property 


roots 


5(x ^^root) 
l/'(a;root)| 


( 102 ) 


one can solve the angular part of the integral if f{x) is taken to be a; ± cuk ~ <^p-k with x = cos 9 
and 6 the angle between the vectors k and p (the -|- sign corresponds to decay and the — sign 
to Landau damping). The sum present in (102) is over the roots of the function f{x), which for 
/(x) = a; ± cuk - Wp^k, are given by 


X = 


=F 2uj Wk 
2|p||k| 


(103) 


with |/^(x)| = |p||k|/c(;p_k. In order to obtain a nonzero contribution, the roots of the function 
f{x) in the d's must be inside the interval [— 1 , 1 ], i.e. 


^ P^ - ± 2a ;ujk ^ ^ 

“ 2|p||k| “ ’ 


(104) 


We analyze the two regions independently: 


1. Decay: In the case of decay the J-function in Eq. (101) implies that a; = a;k -|- Wp-k- This is 
only possible provided to > Vp^ -|- 4m^, so this contribution to damping only occurs above 
the 2-particle threshold. The lower and upper integration limits k- and A:+, which result 
from the restriction (104), can be easily expressed as 


k± = 



(105) 


2. Landau damping: In this case the contribution only occurs below the light-cone (p^ > a;^). 
The integration limits resulting from the restriction (104) turn out to be the same as in the 
case of decay^. Notice that both for decay and Landau damping the function inside the 
square root in the integration limits (105) is positive, so k± are real. 

After the angular integration is performed, the contributions to the imaginary part of the retarded 
self-energy coming from decay and Landau damping can thus be written as 

ImE:^(ca,p) = f {[2nk-M] 0(a;^ - p^ - 4m^)-F 2nk0(p^-a;^)} (106) 

4 Jk- 47r IpI iwk 


This is not true in general. It does not happen, for instance, in the contribution of the sunset diagram to damping 
[ 72 , 78 ]. 


5 















29 


The remaining integrations can be easily performed to obtain 




IGttIpI 


Tin 


1 - e 1 - . ^ 


+ rin 


- 4m^)+ 

( |_e-/3c._ )^(p'-^")V (107) 


with uj± given by 




a; ± IpIW l + 


4j^2 


p2 — 


(108) 


The same result was obtained using Laplace transform methods [79]. We observe from (107) that 
the perturbative retarded self-energy coming from the “eye” diagram does not contribute to on- 
shell damping. The corresponding on-shell plasma excitations (plasmons) are stable and behave 
as free quasiparticles. 

The same conclusion can be obtained by performing the analysis of the damping rate on the 
lattice. In this case, an explicit form for the damping rate such as (107) cannot be given due to, 
among other things, the lack of rotational invariance. The lattice damping rate can be calculated 
by studying the evolution of the mean field, as done in section IV E. 
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